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Abstract 

We present results from a numerical simulation of the two-dimen- 
sional Euclidean Wess-Zumino model. In the continuum the theory 
possesses N = 1 supersymmetry. The lattice model we employ was an- 
alyzed by Golterman and Petcher in where a perturbative proof was 
given that the continuum supersymmetric Ward identities are recov- 
ered without finite tuning in the limit of vanishing lattice spacing. Our 
simulations demonstrate the existence of important non-perturbative 
effects in finite volumes which modify these conclusions. It appears 
that in certain regions of parameter space the vacuum state can con- 
tain solitons corresponding to field configurations which interpolate 
between different classical vacua. In the background of these solitons 
supersymmetry is partially broken and a light fermion mode is ob- 
served. At fixed coupling the critical mass separating phases of broken 
and unbroken supersymmetry appears to be volume dependent. We 
discuss the implications of our results for continuum supersymmetry 
breaking. 
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1 Introduction 



Supersymmetry has often been invoked as a necessary ingredient for any 
particle physics theory which attempts to bridge the gap between the scale 
of electroweak symmetry breaking and the much larger scale associated to 
unification of the low energy gauge interactions. The basic idea is that while 
generic field theories involving scalars are unstable to large radiative correc- 
tions which mix scales, these radiative effects can be made much smaller 
if the scalar theory is embedded inside some supersymmetric theory. The 
dynamical breaking of supersymmetry through non-perturbative effects can 
then occur at scales which are exponentially suppressed relative to the grand 
unified scale. This symmetry breaking can, in turn, then trigger electroweak 
breaking. 

Thus the non-perturbative structure of supersymmetric theories is a sub- 
ject of great interest. The only tool for a systematic investigation of non- 
perturbative effects is the lattice and so a lot of effort has gone into formu- 
lating lattice supersymmetric theories :2 . Typically it is difficult to write 
down lattice actions which can be shown to flow to a supersymmetric fixed 
point, without fine tuning, as the lattice spacing is reduced. 

The model we examine in this paper - the two dimensional Wess-Zumino 
model appears to provide an exception to this rule. This theory involves the 
interactions of scalars and fermion fields and exhibits an N = 1 supersym- 
metry in the continuum. A version of this model defined on complex fields 
and possessing N = 2 supersymmetry was the subject of a recent numerical 
study in j3] and was also examined in a variety of earlier papers [2]. The 
N = 2 model actually possesses an exact lattice supersymmetry which can 
be seen to result from its proximity to a continuum topological field theory 
0. 

We have chosen to study a particular Euclidean lattice formulation of 
the N = 1 model due to Golterman and Petcher pQ. The model has also 
been studied using a Hamiltonian formulation in [H] and [7|. Unlike the 
Hamiltonian formulations, the Euclidean lattice theory does not retain any 
exact supersymmetry. Nevertheless, Golterman et al. prove that the discrete 
analog of the continuum supersymmetric Ward identities are satisfied exactly 
in the limit of vanishing lattice spacing without the necessity of additional 
fine tuning. The proof is perturbative and our goal in these simulations was 
to check whether the model allows for supersymmetry breaking via non- 
perturbative effects. We find that indeed the lattice model shows evidence 
of supersymmetry breaking for small values of the lattice mass parameter. 
Furthermore, this breaking is correlated to the onset of field configurations 
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which sample both the classical vacua of the model. In this limit we also 
observe a light fermion state which we speculate may play the role of a 
Goldstino associated with spontaneous supersymmetry breaking. 

We have developed and tested a fourier accelerated version of the so- 
called R- algorithm [Hj to handle the fermionic integrations. For details of 
this Fourier acceleration technique in the context of the HMC algorithm we 
refer the reader to We have employed an exact algorithm to calculate 
the sign of the Pfaffian resulting from the integration over the fermion fields. 
These issues are discussed in detail in section 2. We present our evidence 
for symmetry breaking together with numerical results on the spectrum and 
Ward identities in section 3. In section 4 we summarize our findings and 
discuss their implications for supersymmetry breaking in the continuum in 
finite and infinite volume. 
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2 Lattice model 



We consider the on-shell two-dimensional Wess-Zumino model represented 
by the following continuum action in Euclidean space [Q: 

s = ( d 2 x I [{d^f + v# + p\m + p 2 m (i) 



where <f> and if) are a real scalar field and a two component Majorana spinor 
respectively. The construction of Euclidean Majorana spinors is described 
by Nicolai in 10 j. The expression Q{4>) = (j) + P'{<p) will be referred to as 
the fermion matrix. The potential P((fi) we consider (actually the derivative 
of the superpotential) takes the following form depending on a mass M and 
a coupling constant G: 



P(<f>) 



M<f> , G = 

G(j) 2 -M 2 /4G , G^O 



Notice that this potential is slightly different from the one considered in 
but may be derived from it by a simple shift in the scalar field. It has the 
advantage that the total action now only depends on M 2 which allows us 
to restrict our simulations to positive M. Notice also that the interacting 
theory has two classical vacua at <p = ±M/2G. The action Q is invariant 
under the following supersymmetry transformation: 

5<f> = Sip = - P{(j)))e 

The simplest super symmetric Ward identity following from this invariance 
takes the form 

< j> s $y > + <[#- P(<p)]x4>y >= (2) 

Integrating out fermion variables in the path integral leads to the follow- 
ing form of the partition function (see the Pfaffian definition in 
appendix) : 

So = /" m GWPfiT^l MHQ T Q)]/±Sb 



Z = J D<f)D^e- So = J D(f> Sign[Pf(CQ)] e 



where C is a Euclidean representation of the charge conjugation matrix and 
Sb stands for the bosonic part of the action: 



S b = J d 2 x \ \(d^f + P' 



4 



In practice we simulate the system without regard to the sign of the Pfaffian 
using the following action S 

S = -hr[MQ T Q)] + jd 2 x l - [{d^f + P\<j>)\ (3) 

The expectation values of physical observables are then obtained by reweight- 
ing with the measured sign of the Pfaffian in the usual manner 

<OSignmCQ^ 
< U >So < Sign[Pf(CQ)] > s { ] 

We now turn to the lattice model. First, we replace the continuum 
derivative operator by the symmetric difference matrix D^ r , where the latter 
is defined as: ^ 

D rr' = ^H-e^r' - <^r-e M ,r'] 

where r, r' are two-dimensional vectors enumerating the lattice sites and 
e„ is a unit vector in the /U-direction jjl = 1,2. In terms of this difference 
operator the fermion matrix on the lattice can be represented as: 

Q = Qrr> = 7^ + 6 ^r' 

where a, f3 are spinor indices. We have employed the following representa- 
tions of the Dirac matrices 

71 = ( -1 ) 72 = ( -1 ) 
The matrix C is given explicitly as 

' -1 \ 



c = 



1 



It is convenient to define an operator n" r ,: 

1 I* 

In particular 

□ 2 r , = (£>"£>")„, 

In terms of the operator n™ r , the lattice potential and its derivative can be 
represented as follows: 



Pr 



m<Pr - □J r /0r'/2 , 5 = 

g<ft - m 2 /4g - , g^O 
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dPy 



mS rr , - Dj r ,/2 
2g<f> r 5 rr , - ul r ,/2 




d(j> T , 



where the term with / is the Wilson mass operator, which serves to 
eliminate problems due to doubling of the lattice fermion modes and vanishes 
in the continuum limit. The dimensionless lattice couplings g and m are 
related to their continuum counterparts through the relations g = Ga and 
m = Ma with a the lattice spacing. 

Finally the lattice representation of the continuum action (j3J) can be 
viewed as the sum of the following boson and fermion components: 



S b = \{-<t> r u 2 rr ,cp Y , + P r P r } 
-itr[ln(Q T Q)] = -i[ln(Q T Q)] 



oca 
rr 
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3 Simulation Details 



To simulate the system (J3J) we use an importance sampling technique based 
on a classical evolution of the fields in some auxiliary time. To implement 
this it is necessary to introduce a Hamiltonian 

H = ^PrPr + S 

associated with this auxiliary time variable t and corresponding momentum 
field p conjugate to the field <j>. On integration over the auxiliary momentum 
p it is trivial to show that the classical partition function associated to 
H reproduces the quantum partition function associated with the original 
action S. The advantage of this Hamiltonian formulation is that it admits a 
classical dynamics, which can be used to generate global moves of the field 
<j>. 

We evolve the system governed by H according to a finite time step 
leapfrog algorithm in the usual manner 

<Pt+st = fa + PtSt + F t (5t) 2 /2 
Pt+st =Pt + (F t + F t+5t )5t/2 

where F is the force associated with the classical evolution. The ergodicity 
of the simulation is provided by periodically drawing new momenta p from 
a Gaussian distribution. In order to decrease the autocorrelation time as- 
sociated with this dynamics we have utilized acceleration techniques similar 
to those explored in j^j- Specifically, the discrete time update of the fields 
corresponding to the Hamiltonian evolution is carried out in momentum 
space with a momentum dependent time step which is tuned so as to evolve 
low momentum components of the field more rapidly than high momentum 
components. Specifically we used a time step of the form 

5t(n)= macc + A 



EjUl sin 2 ^ + (m acc + 2 E'=i sin 



2 

L 



where the lattice momenta n are integer vectors with components ranging 
from — > L — 1 for anlxl lattice. The parameter m acc is typically set to 
the input lattice mass which is close to optimal in these simulations. 

The total force can be represented as a sum of boson and fermion con- 
tributions: 

F r = F r b + F/ 
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The evaluation of the boson force is straightforward: 



F r 6 = -|g = n^-P r ^, r 

In order to evaluate the fermion force we first evaluate the following expres- 
sion involving the fermion matrix: 



db s db r dbs sr ^ sr db r ,db s '" Lyi " c, r 

The fermion force then is: 

P* f = -^ L = \[{Q T Q)-r f {QT d f^ = 9 [{Q T Q)- l rjQ sr 



r's) 



s 



The computation of the fermion force appears to be problematic as it 
appears to require the repeated inversion of the fermion matrix which is 
prohibitively expensive. In order to resolve this problem we use the so-called 
R-algorithm [Hj. The algorithm proceeds by replacing the exact inverse 
matrix (Q T Q)~ 1 by a stochastic estimator given by the following expression 

[{Q T Q)~ l C ^< >n (5) 

where the vector X is defined through a random Gaussian vector R g as: 

QX = Rg 

and the averaging in (JSJ) is accomplished over N different random noise 
vectors R g . 

The larger the number of noise vectors used iV the more accurate is the 
evaluation of the inverted matrix in © , but the longer computational time 
the evaluation takes. It is clear that the optimal value of N is given by that 
which minimizes the error in the inverse matrix for fixed computational time 
T. Defining the norm of a matrix \\A\\ as 



V l 3 

The relative error is then 

5\\A\\ IN (5\\A\\ 



\A\\ V T I IL4 



N 
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where {5||^4||/||A||}jv is the relative error produced by a single application 
of an R-algorithm with averaging over ./V noise vectors. Hence the relative 
error obtained over time T can be characterized by the algorithm efficiency 
E which we define as 

■S\\A\ 



E = y/WU 



\A\ 



N 

Our tests showed that this algorithm efficiency does not depend strongly on 
the choice of N (Figure Furthermore, we monitored the average bosonic 
action < Sb > and observed no systematic drift with N. Consequently we 
chose N = 1 in all our runs. In this limit the corresponding fermion force 
term yields 

F f w i x « x e W T d Q)R = lqg(x«Xg + xfx?) 

Finally let us turn to the issue of the sign of the Pfaffian which results 
from integrating out the Majorana fields. As we have stressed the simula- 
tion action discussed above utilizes the absolute value of this Pfaffian and 
observables must be re- weighted by the sign of the Pfaffian in order to com- 
pute physical expectation values. We have chosen to use an exact algorithm 
to compute this sign. Since we are in two dimensions and need only do 
this when making measurements this turns out to be quite manageable in 
a practical sense. Our procedure was outlined in and for completeness 
we list the proof and details of the algorithm in an appendix. In essence the 
original antisymmetric matrix can be transformed to a special block diago- 
nal form via a similarity transformation built from a triangular matrix. The 
determinant of the latter can be shown to yield the Pfaffian. We then fold 
the sign of the Pfaffian in with measurements of observables according to 
©. This reweighting procedure is an effective way to measure expectation 
values of a variety of observables. However, in certain conditions this tech- 
nique may fail. The following arguments highlight the problems that may 
be encountered in this type of situation. 

Let N + and N_ be the numbers of configurations with positive and 
negative values of Pf(CQ) obtained from the simulation of the system (j3J). 
Then the average value < O >s of any physical observable O in the system 
Q can be evaluated using @ as: 

where 0± are average values of O obtained in configurations subsets with 
positive and negative Pf(CQ). This averaging procedure reveals two statis- 
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tical problems. The first problem is that if iV + ~ iV_ (that is the probabil- 
ities to find the system with either sign of the Pfaffian are approximately 
the same) then the error of the evaluation ((HJ) experiences an amplification 
by a large factor (N + + N_)/(N + — 2V_). In this case it is possible for the 
error to swamp the signal in the measured value < O >. Although acquir- 
ing more measurements will decrease the fluctuations it might not solve the 
amplification problem if 

Jim —t~l (7) 

A second problem is that the expression © provides a good evaluation 
of < O >s Q only if < O > is uniquely defined in the limit N + + iV_ — ► oo, 
which is not necessarily the case. If (JJJ) takes place then < O > is well 
defined only if < + >=< 0_ >, which is not guaranteed to be true. 

In practice we find that many of our observables suffer large and difficult 
to quantify errors for small values of the lattice mass where we observe oscil- 
lations in the sign of the Pfaffian. This precludes making strong quantitative 
statements in that region. 
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4 Results 



We obtained data for lattice sizes L = 8 and L = 16 for a fixed lattice cou- 
pling g = 0.125 while varying the lattice bare mass m. The classical vacua 
of the lattice theory correspond to vanishing fermion field and boson field 
(ft = ±m/2g. For large m, field configurations which interpolate between 
these two vacua are associated with large values of the action and are hence 
expected to be highly suppressed. We thus expect the boson field to be con- 
fined in the neighborhood of one of the classical vacua for sufficiently large 
mass. In the continuum the action is invariant under (ft — > —(ft implying that 
these two vacuum states are equivalent. This is no longer true on the lattice 
due to the presence of the Wilson term (actually the sign of the Pfaffian 
may also change under this symmetry). Indeed our simulations reveal that 
only the state with < (ft >~ —m/2g survives at large m. As m decreases 
we expect that tunneling to the other vacuum state may occur and this is 
indeed seen in our simulations. Figures |2] and 01 show plots of < Pf(CQ) > 
and < (ft > /m versus m for L = 8, 16. Below some critical m = m c (L) the 
sign of the Pfaffian which is negative at large mass m starts to fluctuate. 
Additionally, in this region we can see that the average field < (ft > /m also 
undergoes large fluctuations which are the direct result of the Pfaffian sign 
changes. Indeed, at small mass we observe that for each configuration in 
our ensemble the sign of the Pfaffian is very accurately correlated with the 
sign of the mean boson field. Figure 0] shows a time series of both quantities 
at m = 0.125 and L = 8 which illustrates very well this behavior. Actually 
it is easy to see why this is so. Imagine expanding the Pfaffian as a power 
series in the boson field (ft. For sufficiently small m we expect that only the 
leading term is important and by translational symmetry this can depend 
only on the field summed over all lattice sites. 



As we discussed in the previous section this sign oscillation renders accurate 
measurements of < (ft > and its error very difficult in this region. 

We have also measured the (zero momentum) boson and fermion cor- 
relation functions over the same range of lattice bare masses. Figures 
and show typical bosonic and fermionic two point functions computed on 
ensembles corresponding to L = 16 with m = 0.5. These are fitted by hy- 
perbolic cosh and mixed hyperbolic sinh and cosh functions to extract the 
corresponding boson and fermion masses. These (lattice) masses are shown 
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N + 


N_ 


40968 


59032 


43814 


56186 


52252 


47748 



Table 1: Numbers of positive and negative Pfaffians for L = 8 and m = 0.25 



in figures [7| and |H] for L = 8 and L = 16 respectively. The statistical errors 
we show neglect the effects of correlation between observables at different 
timeslices. Consider the data for L = 8. Notice that the boson and fermion 
masses are equal within statistical errors at large bare input mass but de- 
viate substantially at small mass - the diagonal spinor components of the 
fermion correlator being dominated by a light state. Contrast this will the 
off-diagonal components of the fermion correlator for small bare mass which 
yield a much heavier mass degenerate with the boson mass within statistical 
errors. A light fermion state is also visible in the L = 16 data at small mass. 
The mass of this light fermionic state appears to decrease with the bare 
input lattice mass m. It is tempting to conclude from these observations 
that for small enough mass supersymmetry breaks as a result of mixing be- 
tween the two classical vacua - this being signaled by the appearance of a 
Goldstino. 

Another line of evidence in favor of this derives from the partition func- 
tion itself. On a finite lattice equipped with periodic boundary conditions, 
such as employed in our simulations, the partition function can be thought 
of as yielding a representation of the Witten index. Vanishing Witten index 
is a necessary condition for supersymmetry breaking. But Z can also be 
related to the expectation value of the sign of the Pfaffian in our simulation 
ensemble 

Z So =< Sign (Pf(CQ)) > s 

Thus we see that a vanishing partition function would require equal numbers 
of positive and negative sign Pfaffians in our ensemble. Table 1. shows the 
numbers of positive N + and negative iV_ Pfaffians for three different runs 
at the same parameter values L = 8 and m = 0.25 each containing 100, 000 
measurements. While the relative errors of on the order of ten percent it 
should be clear that the data are consistent with a vanishing Witten index 
implying a non-zero vacuum energy. 

To investigate this symmetry breaking further we have looked at the 
simplest super symmetric Ward identity involving two point functions ((2J). 
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Figures to ^] show the bosonic and fermionic diagonal and off-diagonal 
spinor contributions to this Ward identity together with their sum for two 
different values of the bare mass m = 0.125 and m = 0.5 on a lattice with 
L = 16. Clearly for large mass this relation is satisfied within errors for 
all spinor channels but it is clearly violated at small mass for the channel 
involving the diagonal spinor correlations. The latter channel is precisely the 
one in which the light fermion was seen and support the idea that breaking 
of supersymmetry is associated with the appearance of a Goldstino. Notice 
that the off-diagonal components of the Ward identity are still accurately 
satisfied even at small mass. We will argue in the next section that this 
is exactly what we might expect for a partial breaking of supersymmetry 
associated with the appearance of a finite volume vacuum state composed 
of solitons. 
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5 Conclusions 



We have studied a lattice regularized version of the two-dimensional Wess 
Zumino model which possesses N = 1 supersymmetry in the naive contin- 
uum limit. This model was first analyzed in [I] where it was shown per- 
turbatively that the supersymmetric Ward identities are recovered without 
finite tuning in the limit of vanishing lattice spacing. The goal of our simu- 
lations was to check these conclusions at the non-perturbative level and to 
specifically to address the important issue of supersymmetry breaking. We 
have considered the model for fixed lattice coupling g = 0.125 and varying 
lattice mass m for two lattice sizes L = 8 and L = 16. For large m our 
results favor a supersymmetric phase in which boson states pair with equal 
mass fermion states and the supersymmetric Ward identities are satisfied. 
In this region of parameter space corresponding to <j) ~ —m/2g the boson 
field suffers small fluctuations around a single vacuum state. 

As the mass is lowered however this picture changes and below some crit- 
ical mass m c (L) we see configurations in which the mean boson field varies 
in sign corresponding to tunneling between different vacua in auxiliary time. 
The appearance of states which interpolate between different perturbative 
vacua is of course entirely a non-perturbative effect. Associated with these 
tunneling states we see oscillations in the Pfaffian of the fermion operator 
and the appearance of a light fermion visible in the diagonal components of 
the fermion correlator. In such a phase it appears that supersymmetry is at 
least partially broken. 

It is possible to get some further understanding of this phenomenon 
within the context of the semi-classical approximation. Consider first the 
continuum model. It is clear that in a finite volume corresponding to a box 
of size Lphys) i n addition to the supersymmetric vacua <j) = ±M/2G, there 
are additional local minima of the action (^Q) corresponding to domain wall 
solutions which interpolate between these vacua. 

4>{x) — > A, x — > oo 
4>(x) — > —A, x — ► — oo 

where A = M/2G and x corresponds to one of the coordinate directions. 
Indeed, in the continuum, these solutions take the form 

<j){x) = Atanh(M(x - x )/2) 

While the mass of such a soliton state is non-zero it is possible to show 
that it is nevertheless annihilated by a single component of the Majorana 
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supercharge and hence such a state preserves one half of the original super- 
symmetry ^3] . This is the origin of our observation that certain components 
of the Ward identities appear to be satisfied at all values of the parameters. 
The action of such a soliton is easily evaluated «Sbw = zGL p h ys A 3 and be- 
ing proportional to the integral of a total derivative term is topological in 
character. The corresponding free energy associated with such domain wall 
solutions then varies as 

Fdw ~ — InLphys + Sow 

where the logarithmic variation with volume arises from the number of ways 
the domain wall can be introduced into the finite volume. These arguments 
lead one to conclude that these non-supersymmetric vacua will dominate 
over the supersymmetric vacua if 

■Sew = ( hi L phys 

Lphys 6G 2 I -^phys 

At fixed G this result is in qualitative agreement with our lattice results since 
it predicts a critical mass Mc(-L p h ys ) below which supersymmetry would be 
broken. Translating this result naively into lattice variables leads to the 
prediction that mc ~ 0.3 for L = 8 and g = 0.125 which is quite close to the 
continuum estimate Mca = 0.46 for L p h ys = 8a. According to our observa- 
tions this critical mass shifts to smaller values as the lattice size increases 
which is also in agreement with these analytic arguments. Furthermore, in 
the vicinity of such a domain wall the fermion is approximately massless and 
so can play the role of a Goldstino associated with supersymmetry breaking. 
Notice that these arguments rely on the constraint of finite volume. The ac- 
tion of such a soliton is unbounded in infinite volume and hence we would 
naively expect solitons to be completely suppressed in such a limit. 

Of course we would like to know whether this finite volume supersym- 
metry breaking scenario persists in the continuum limit. In general, in finite 
physical volume V = L^ hy , the continuum limit a = 1/L p h y — ► should be 
approached by fixing (in this case) two renormalized physical parameters 
which may be taken as the mass MnL p \ iy and coupling constant GrL^ - 
expressed in units of the physical length L p h y - Perturbation theory allows 
us to relate these renormalized dimensionless quantities to their bare lattice 
counterparts 

G R L phy ~ gL, M|L^ hy ~ m 2 L 2 - Cg 2 L 2 In (jm) 
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where fi is the mass scale associated to the renormalization point and C is 
a numerical constant. Along such an RG trajectory the value A = ^ can 
then be related to its constant (continuum) value A^j via the relation 

^(i-Hi))* 

For small enough L < L c the log term on the right is small and A ~ A^j. If 
this is the case the soliton action S*dw ~ 9^X 3 is approximately constant and 
the corresponding free energy of such configurations negative for any value of 
the bare parameters. This would correspond to finite volume supersymmetry 
breaking. Conversely for large enough L the log term will dominate and lead 
to an infinite soliton action as L — > oo for any value of the bare parameters. 
In this limit the solitons should disappear and supersymmetry should be 
restored. The cross-over between these two behaviors occurs for 

L ~ L c = Ae~ x2 

These conclusions are in agreement with the reasoning presented in |14j . 

While this work was in preparation we received a preprint ^H] m which 
the same model is studied in a Hamiltonian framework. The conclusions of 
this study are broadly in agreement with ours. 
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Appendix: the algorithm for determining the Pfaf- 
fian of an antisymmetric matrix. 

In this chapter we describe the algorithm for determining the Pfafhan of an 
arbitrary antisymmetric 2N by 2iV matrix M, which is defined as follows: 

PfM = £a iA,-,aAr,/3Ar M ai,/3i ■ ■ ■ M a N ,/3 N (9) 

The algorithm utilizes the following theorem. 

Theorem. If P is a matrix such that an antisymmetric matrix M can be rep- 
resented as M = P T JP where J = diag(z73, 273, 473) is a block-diagonal 
matrix then PfM = det P (here C = 273 is the Euclidean representation of 
the charge conjugation matrix for two-dimensional system). 

The theorem can be proved using the representation of Pfafhan in terms 
of an integral over a Grassmann 2N-vector 9. Defining 9 = P9 we have: 

PfM = J d9e~¥ T Me = J dde -^p^jpe = J d6e -\e T Je = 

= J d9— [9~2n-l92n] N = J d9— [P2n-l,a^aP2n,/3&p] N = 

= J d9Pi tai 9 ai P 2) p 1 9[ 3l . . . P2N-l,a N 9a N P2N,f3 N 9/3 N = 
= £a 1 ,/3 1 ,...,a N ,l3 N Pl,a 1 P2,l3 1 ■ ■ ■ ?2N -l,a N P2N ,f3 N = det P 

Notice that the matrix P is not orthogonal (P T 7^ P^ 1 ), hence it is not 
associated with any basis transformation in 2iV-dimensional vector space. 

The above theorem can be given an alternative formulation. Defining 
Q = P' 1 leads to the following statement: if Q T MQ = J then PfM = 
(detQ) -1 . This formulation is used in the algorithm we describe below. 
The purpose of the algorithm is to represent a given antisymmetric matrix 
M in terms of a triangular matrix Q so that the det Q and hence the PfM 
can be found easily. 

The algorithm task. Given an arbitrary antisymmetric 2N by 2N matrix M 
find a triangular matrix Q such that Q T MQ = J, J = diag(i^,i^, 273). 
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The triangular matrix Q = {qi} represented in terms of its columns % 
will satisfy the relation above iff its columns satisfy the following conditions: 

(ga-lM^-i) = 0, (q 2i Mq 2j ) = 

(Q2iMq 2j -i) = ~{q 2 j-iMq 2i ) = <% 
The following algorithm by construction leads to such a matrix Q. 

The algorithm. 

1. Establish a unary 2N by 2N matrix Q = diagjl, 1, 1} = {ej}, where 
unary vectors are columns of Q, i = 1,2, ...,2N. 

2. For odd values i = 1,3, 2N — 1 repeat the following steps: 
2-a. Leave ej as is. 

2-b. Redefine e*+i — ► ei+i/(ej+iMej) 

2-c. For k = i+2, i+3, 2N redefine — > e^ — ej(ei + iMefc)+ei + i(eiMeA;) 
Notice that in 2-c the vector ej+i is used after it is redefined in 2-b. 
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Figure 1: Algorithm efficiency as a function of the number of noise vectors. 
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Figure 2: Average field and average Pf(CQ) for L=8. The field values are 
rescaled for m < 0.46 by a factor of 1/200. 
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Figure 3: Average field and average Pf(CQ) for L=16. The field values are 
rescaled for m < 0.29 by a factor of 1/8. 
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Figure 4: Evolution of Pf(CQ) and < > in auxiliary time t for L=8, 
m=0.125. 
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Figure 7: Mass gaps for L=8. Mass gaps from bosonic correlators are shown 
for m < 0.46 without error bars 
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Figure 8: Mass gaps for L=16. 
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Figure 10: Contributions to the off-diagonal components of Ward identity 
for m=0.125. 
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Figure 11: Contributions to the diagonal components of Ward identity for 
m=0.5. 
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Figure 12: Contributions to the off-diagonal components of Ward identity 
for m=0.5. 
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